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We propose to combine the Floquet formalism for systems in ac fields with the dynamical mean- 
field theory to study correlated electron systems periodically driven out of equilibrium by external 
fields such as intense laser light. This approach has a virtue that we can nonperturbatively include 
both the correlation effects and nonlinear effects due to the driving field, which is imperative in 
analyzing recent experiments for photoinduced phase transitions. In solving the problem, we exploit 
a general theorem that the Hamiltonian in a Floquet matrix form can be exactly diagonalized 
for single-band noninteracting systems. As a demonstration, we have applied the method to the 
Falicov-Kimball model in intense ac fields to calculate the spectral function. The result shows that 
photoinduced midgap states emerge from strong ac fields, triggering an insulator-metal transition. 
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I. INTRODUCTION 

Controlling phases of matter is a central issue in the 
physics of strongly correlated electron systems, where a 
rich variety of phases are realized for various physical 
degrees of freedom such as spin and charge. Since a 
phase transition dramatically alters macroscopic proper- 
ties of the system, it is of great importance to know how 
the phase changes as external parameters (temperature, 
pressure, band filling of the system, etc.) are varied, in 
equilibrium. 

Now, the last decade witnessed that a controlling 
factor, intense laser fields, can trigger "phase transi- 
tions" in correlated electron materialsf^i^ One repre- 
sentative class of materials is perovskite manganites, 
in which an insulator-to-metal transition is induced by 
photoexcitationi^i^ Recent experiments also indicate that 
a ferromagnetic spin alignment can emerge in the in- 
duced metallic phase in manganites Such phenom- 
ena are called photoinduced phase transitions (PIPTs), 
where irradiation of photons allows one to change elec- 
tronic, magnetic, optical, or structural properties of the 
system. 

The photoinduced phase transition, however, is dis- 
tinct from conventional phase transitions in equilibrium 
in that the photon field drives the system out of equilib- 
rium. In other words, the change in phases in nonequi- 
librium challenges our understanding of phase transitions 
which is normally conceived in equilibrium. In the PIPT, 
we have to consider, on top of the nonlinear electric-field 
effect, the electron correlation effect. The nonlinear effect 
appears as a threshold behavior, i.e., a macroscopic tran- 
sition only occurs when the intensity of the driving field 
exceeds a certain strength. On the other hand, important 
correlation effects such as Mott's metal- insulator transi- 
tion are nonperturbative effects. Thus, we cannot employ 
the linear-response theory^ nor a mean-field treatment of 
the electron-electron interaction. 

Here we propose a theoretical approach^ for photoin- 
duced phenomena, where the Floquet-Green function 



method (FGFM)^i^ d°d-'^ is plugged into the dynamical 
mean- field theory (DMFT)ii2i This formulation provides 
us with a promising way to treat both of the nonlinear 
effect of the electric field and the correlation effect simul- 
taneously in a nonperturbative manner. We then apply 
the method to study the response of the Falicov-Kimball 
(FK) model, one of the simplest lattice models for 
correlated electrons, to ac electric fields. A particular 
emphasis is put on the technical basis of FGFM. 

During the preparation of the present study, we notice 
that Joura et ali^ adopted a technique similar to FGFM 
to solve the Dyson equation for a system in dc fields. 
Here we present a general framework of DMFT out of 
equilibrium for arbitrary time-periodic fields in a more 
transparent viewpoint exploiting the Floquet formalism. 
In this context we should emphasize that FGFM is not 
just a numerical technique but offers a fruitful physical 
picture for nonequilibrium systems as revealed here. 

The paper is organized as follows: in Sec.|TTl we review 
the Floquet theorem and the Floquet matrix method, on 
which our theoretical description is based. The rest of the 
paper is devoted to our original formulation. In Sec. IIIII 
we define the Floquet matrix form of Green's function, 
which is our starting point of FGFM. Then in Sec. IIVI 
we derive a general expression for Green's function and 
its inverse for noninteracting electrons in a Floquet ma- 
trix form. We calculate noninteracting Green's functions 
for several cases to discuss their physical implications. 
In Sec. lYl we derive a general theorem that identifies the 
eigenvalues and the eigenvectors of a Floquet matrix form 
of the Hamiltonian for noninteracting electrons. While 
Sees. lIVl and IVl address noninteracting cases, we move on 
to correlated electron systems by incorporating FGFM 
in DMFT in Sec.|Vll We formulate the obtained Green's 
function in a gauge-invariant manner in Sec. IVIII We 
then apply our method to the FK model in Sec. IVIIIl 
and calculate the spectral functions for dc and ac fields, 
where we discuss how the Mott-insulating state is trans- 
formed into the metallic states in the external fields. We 
summarize the paper and give future problems in Sec. lIXI 
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II. FLOQUET THEOREM AND FLOQUET 
MATRIX 

An external field drives an electron having an energy e 
into another state with a different energy, where there are 
many scattering channels in nonequilibrium. However, if 
the driven system is periodic in time with a frequency 
f2, the allowed channels are limited to the processes such 
that e e + nQ, where n is an integer. This greatly 
reduces the channels' degrees of freedom to be dealt with. 
One can take advantage of such a simplification through 
the Floquet matrix methodji^iiiii^ of which we give an 
overview in this section. 

The method has been used as an effective approach 
toward photoexcited systems. The concept of the Flo- 
quet matrix originates from the Floquet theoreroi^ for a 
periodically driven system, an analog of the Bloch theo- 
rem applied to a spatially periodic system. Floquet the- 
orem is a general theorem for differential equations of a 
form dx{t)/dt = C(t)x{t) with C periodic in i, which 
include equations of motion for systems subject to exter- 
nal driving forces that periodically oscillate in time. One 
representative example is the Mathieu equation which 
describes a parametric resonance phenomenon. Here we 
restrict ourselves to a quantum system whose dynamics is 
determined by the time-dependent Schrodinger equation, 

zj^^it) = nmit), (1) 

where is a state vector of the system, and H{t) is 
the time-dependent Hamiltonian, which is assumed to be 
periodic in t, 

H{t + T)=H{t), (2) 

with a period r. The Floquet theorem states that there 
exists a solution of Eq. ([T]) which is an eigenstate of the 
time translation operation t —^ t + t, implying 

^^,^{t) = e-*^°*u„(i) (3) 

with e"'^"'^ as an eigenvalue of the time translation, a 
as a set of quantum numbers, and Ua(t) = Ua{t + r) as a 
periodic function of t. Hence we can Fourier expand Ua {t) 
as Ua{t) = X]„e~™^*u2 with the frequency Q — 2t:/t, 
where u" is called the nth Floquet mode of Floquet state 
([3]). We can then Fourier transform Eq. ([T]) to have 

^-ff™n< = (ea+mr^X, (4) 

n 

where 

Hrnn = - f dt e*('"-")"*ij (t) (5) 
T J-t/2 

is the Floquet matrix form of the Hamiltonian. The fac- 
tor Ea + TTi^ appearing on the right-hand side (rhs) of 
Eq. ^ is called quasienergy, which forms a ladder of en- 
ergies with a spacing 17. Since the Hamiltonian is time 



dependent, the energy is not conserved in general. How- 
ever, Eq. (dl shows that the energy is conserved up to an 
integer multiple of fi, corresponding to the absorption or 
emission of the photon with the energy SI. Each element 
in the Floquet matrix Hmn corresponds to the proba- 
bility amplitude of the transition from the mth Floquet 
mode to the nth one, so that off-diagonal components 
represent excitations driven by the external field while 
the diagonal ones the probability to remain in the same 
mode. 

The consequence of the Floquet theorem is remark- 
able: Eq. (dl resembles the static Schrodinger equation in 
equilibrium except for the presence of the Floquet mode 
index n, which means that we have no longer to solve 
the time-dependent Schrodinger Eq. ([T]), in favor of the 
time-independent Eq. ^ . This is the great advantage of 
the Floquet matrix method, which also plays a crucial 
role in Grecn's-fmiction approach. 

III. FLOQUET REPRESENTATION OF 
GREEN'S FUNCTION 

Besides the Floquet analysis of the Hamiltonian, 
we can alternatively describe nonequilibrium states 
in Green's-function approach based on the Keldysh 
formalism. ■^"'^^ The approach of the Floquet matrix 
proves its own worth when it is used within Green's- 
function formalism. The idea of FGFM was first intro- 
duced by Faisal,^ followed by several groups i^iiSiiii In this 
section we give another way to define a Floquet matrix 
form of Green's function, which we shall use in this paper. 

A Green's function has two independent arguments of 
time, t and t' , as denoted by G{t, t'). We define variables 
trci = t — t' and tav = (t+t')/2. In equilibrium the system 
is invariant against continuous time translation, so that 
Green's functions depend on t and f only through iid, 
which enables us to Fourier transform them into func- 
tions of Lu. However, when the system is driven out of 
equilibrium, they generally depend on both tj-ci and fav 
Since the periodic system that we consider in this pa- 
per has the discrete time translation invariance [Eq. ([2])], 
it is guaranteed that Green's function is also invariant 
against tav t^^v + t. For an arbitrary function G{t,f) 
(not limited to Green's function) that satisfies the peri- 
odicity condition, G{t + T, t' + t) — G{t, t'), we can define 
the Wigner transformation of G as 

J-oo T J-t/2 

(6) 

We call G„(w) the Wigner representation of the func- 
tion G. Using the Wigner representation, we define the 
Floquet matrix form of G as 

Gm„(w) = Gto_„ ( w H — ill, (7) 

and call Gmn{^) the Floquet representation. Hereafter 
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a function with one index n should be understood as a 
Wigner representation, while two indices m and n mean 
a Floquet representation. In the Floquet representation, 
we use the reduced zone scheme, i.e., the range of lo is 
restricted to the "Brillouin zone" on the frequency axis: 
—£7/2 < oj < il/2. We can readily check that definition 
is equivalent to the one given by Refs. [sllgl flollTll The 
Floquet representation is used during calculations, while 
the Wigner representation is used when we interpret the 
result, since the Wigner representation has a clear phys- 
ical interpretation that G„ is the nth oscillating mode in 

tav of G(t,i')- 

Actually, the Floquet representation Gmn{'^) has a 
one-to-one correspondence with the Wigner representa- 
tion Gi{uj'). Gi{uj') — > Gjnni^)'- the integers m and n 
should obey the conditions 



n 



n = e 

m 



2 - 2' 



(8) 
(9) 



for H. and lo' . There are two consecutive integers k and 
k + 1 which can be equal to m -I- n satisfying Eq. 
Either fc or A: -I- 1 is congruent to £ modulo 2. Thus m -t- 
n is uniquely determined via Eq. ([9]) and the condition 
m + n = m — n = £ (mod 2). Together with Eq. ([8]), 
TO and n are uniquely determined. For such to and n, 
Lu is given by lu' — (to 4- n)fl/2. Gmn{^) ^ Gi(uj'): t 
and w' are uniquely determined by = to — n and lo' = 
w (to n)f7/2. 

We can immediately realize the advantage of the Flo- 
quet representation in multiplications of two Floquet- 
represented functions. As shown in Appendix |X1 the 
mapping from G{t,t') to Gmn(^) preserves a multipli- 
cation structure, 



dt" A{t,t")B{t" ,t') = C(t,i'), 
^^Ami{uj)Bin{u}) = Gm„(w). 



to be taken in numerical calculations is usually small be- 
cause sufficiently high-order processes should tend to be 
irrelevant when the driving field is not so large, which 
also supports the usefulness of FGFM. 



IV. NONINTERACTING ELECTRONS 

Having defined the Floquet representation of Green's 
function in Eq. ([Tj), we then compute the Floquet- 
represented Green's function for noninteracting electrons. 
Although FGFM has been used by several authors to 
study noninteracting electrons driven out of equilibrium, 
there are still further developments yet to be explored. 
This has motivated us to present an exact and unified 
description of FGFM that can be applied to general non- 
interacting single-band systems in this section. 

In Sec. IIV Al we provide a general expression of the 
Floquet representation of Green's function for any single- 
band model. Then in Sec. IIV 51 we derive the inverse of 
Green's function, which will be used to build DMFT in 
the Floquet matrix form in Sec. IVII After that, we show 
several examples of Green's function for the hypercubic 
(Sec. W^ and other lattices (Sec. HVDj) . Through- 
out the paper we restrict our discussion to a single-band 
model, and omit spin degrees of freedom for simplicity. 



A. General lattices and fields 

Let Cfc be a band dispersion of the system. We make 
the system subject to a homogeneous time-dependent 
electric field periodic in t. Here we choose the temporal 
gauge or the Hamiltonian gauge in which the scalar po- 
tential = 0. Replacing the momentum k with k~eA{t) 
[A{t): a vector potential] in gives the noninteracting 
Hamiltonian, 



H{t) = ^efc-eA(t)4cfc: 



(11) 



As an example, the Floquet representation of the Dyson 
equation is simplified into 

(Gfc)m„(cj) = (G^)m„(w) -I- ^ (G^)mm'(t^)(Sfe)m'„'(cj) 

m'n' 

X {Gk)n'n{0j), (10) 

where G and G° are, respectively, the full and the non- 
interacting Green's functions and E is the self-energy. 
Note that each function has the additional 2x2 matrix 

structure, G = ( „ ) in the Keldysh space (with 



G' ^ 

the three linearly independent components: the retarded, 
the advanced, and the Keldysh one). Thanks to the usual 
multiplication rule of the linear algebra, one can solve the 
Dyson Eq. ^ as Gk{uj) = [G^"\w) - Efc(w)]-i. In ad- 
dition, a typical size of a Floquet matrix that is needed 



where c\, and Ck are the creation and the annihilation 
operators of the electrons, respectively, and we treat the 
electric field as a classical one. The retarded Green's 
function for noninteracting electrons reads 

G^\t,t')^-ie{t-t'){[cu{t), 4(t')]+)o 

= ~ie{t ~ t') exp ^-i dt" [efe_eA(t") - m] 

(12) 

where 9{t) represents the step function, [ , ]+ is the an- 
ticommutation relation, (•••)o is the statistical average 
with respect to the initial density matrix po — e^'^^'^^^"-' 
(where the system is assumed to be in equilibrium with 
the temperature j3~^ at t — — oo), and ^ is the chemical 
potential of the system. We can transform Eq. into 
the Wigner representation via Eq. ^ , and then into the 
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Floquet representation through Eq. The details of 
the calculation are described in Appendix |B] The final 
result is 



(Gr)™nM=E 



1 



uj + £Q + (efc)o + it] 



dx 



27r 



X exp 



dz [ek-eA(z/a) - (efe)o] 



(13) 



where 77 is a positive infinitesimal, and (efc)o is the time- 
averaged dispersion, which coincides with the zeroth Flo- 
quet mode of e^. 



}m — n 



dz 



i(m—n)z /-I ^\ 

^ e ^ ' ek-eA{z/n)- (14) 



Equation (|T3|) is the general Floquet representation of 
Green's function for the noninteracting system driven by 
a periodic field. What is notable about expression (fT3)) 
is that it can be decomposed into well-behaved matrices. 
Let us define two Floquet matrices, 



(Afc), 



dx 
2^ 



exp (^--^ J [<^k-eA(z/n) ~ (efe)o] 



and 



{Qk)mn{^) — 



+ nil + /i - (efc)o + irj 



(15) 



(16) 



where is unitary as shown in Appendix [Cl and Qki^^) 
is a diagonal matrix. The physical meaning of these ma- 
trices will be given later in Sec. |Vl Here let us just note 
a strikingly simple decomposition, 



G^''iLo)^Ak-Qk{Lo)-Al 



(17) 



where we denote a multiplication of a Floquet matrix by 
" •," and omit Floquet indices in Eq. ([T7]). The decompo- 
sition [Eq. (fT7|) ] is essentially used to derive the inverse 
of Green's function in Sec. IIVBI 

The Floquet representation of the advanced Green's 
function is equal to the Hermitian adjoint of the retarded 
one: (Gf )„„(t^) = (Gf^)„„(u;) = (G«");„(a;). Using 
Eq. (HZl), we have G^O(w) = A^ • qI{uj) ■ A^. 



B. Inverse of Green's function 

When one solves a Dyson equation such as Eq. (fTO]) 
to include effects of interaction, the noninteracting part 
appears as an mverse, G^° (uj), rather than Gf"(w) 
itself. Using relation (fT7|) and the unitarity of A^ , we can 



analytically invert Green's function as G™ ^ {us) = A^ • 
Q'i^^(uj) • A^, which can be evaluated exactly as presented 
in Appendix [D] The derived expression for the inverse of 
Green's function is 

(Gfc° ^)„i„(cj) = (cc) + 7iri + ^ + i77)^„i„ - (efe)„i_„. (18) 

One can also prove Eq. from the Dyson equation 
for G™(t, t') in a straightforward manner. Relation (fTS]) 
means that Green's function is the kernel of Eq. ([4]), 
or that the Floquet representation of Green's function is 
equivalent to the inverse of the Floquet matrix form of 
the quasienergy minus the Hamiltonian. One can use 
Eq. (|18[) for any single band Hamiltonian with a homo- 
geneous electric field periodic in time. In Sees. IIV Gl and 
IIVDI we present examples of the calculations that utilize 
relation P^ . 



C. Hypercubic lattice 

As a first example, let us consider a simple cubic lattice 
in d dimensions, whose energy dispersion is given by 



-2t cos ki 

i=l 



(19) 



where t is the hopping and we set the lattice constant 
a — 1. For simplicity we assume that the vector poten- 
tial A{t) is parallel to (1, 1, . . . , 1) with each component 
Ai{t) = A{t). Substituting with ki-eA{t) in Eq. p^ . 
we have 

4-eA(t) = COS eA{t) -f sin eA{t) , (20) 
where we have defined 



'2t sin ki 



(21) 



after Turkowski and Freericks;^ Note that e^^ has the 
odd time-reversal symmetry. Every equation including 
e^j^ and e^'^ must be consistent against the time-reversal 
operation. For instance, one usually finds the factor e^^ + 
ie|f , which is time-reversal even, since the imaginary unit 
i is time-reversal odd. 

An integral over k is performed through 



(22) 



which is called the joint density of states (JDOS).^^ This 
contrasts with the equilibrium cases in which an inte- 
grand depends on k only through el*^, so that we can 
replace the integral variable k with e^^ accompanied by 
the usual density of states p(e). The analytic expression 
for JDOS in arbitrary d dimensions is summarized in Ap- 
pendix[El In particular, in infinite dimensions the JDOS 
becomes a Gaussian function [Eq. (|E2p ].^^ In the follow- 
ing, we consider two kinds of electric fields, a dc field 
rSec. lIVGll) . and an ac field f Sec. IIV G 21) . 
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1. Hypercubic lattice in a dc field 

A homogeneous dc field is given by the vector potential 
proportional to time, 

eA{t)a = nt, (23) 

where 

n = -eEa (24) 

is the Bloch frequency. The system with the field 
[Eq. ((231)] fulfills the periodicity condition [Eq. ^] with 
the period 27r/fJ because of the periodic potential of the 
lattice. With Eq. (fT3|) . Green's function has the following 
Floquet representation, 

(0^.(.)^e-(—)^-E . + + , + 

xJ_.(|)j.-„(|), (25) 
where J„(z) is the nth-order Bessel function, and 

a - \/(^)^+^' (26) 

tan0fe = el7e|^ (27) 

To see how Green's function ((25l) behaves, we calculate 
the local spectral function A„(tj) = ^ ^Im ^j^(G'™)„(a;) 



in infinite dimensions. First, we transform Eq. (|25p into 
the Wigner representation, and then integrate it over k 
with the JDOS [Eq. (jE2|l ]. After taking the imaginary 
part, we arrive at 

A„{uj)^d„^o^aiS{uj + m + fi). (28) 

e 

Here the coefficients in front of the delta functions are 
ai — e~^/^^ /f(l/2il^), where Ii{z) is the modified 
Bessel function of the first kind. Note that the coeffi- 
cients satisfy the normalization condition, ai = 1. 

Equation ([^5)) shows that the local spectral function 
has only the n = component, indicating that the spec- 
tral function evolves into a time-independent function 
[Eq. for a sufficiently long time elapsed after the 

dc field began to drive the system32> Later in Sec. IVII AI 
we shall show that the components with n vanish 
due to a symmetry in the system. 

The spectral function on the hypercubic lattice, dis- 
played in Fig. m consists of a set of delta functions with 
a spacing fl, namely, the Wannier-Stark ladder. The 
width of each peak (approximately effective hopping) is 
infinitesimal due to the Bloch oscillation, where an elec- 
tron is not free to run in a lattice, but only oscillates with 
the frequency Q. 

The inverse of the Floquet-represented Green's func- 
tion in this case is given via relation ((TH]) by 



[Gk )mri(t^) (t^ + nCl + /i + ir])6„ 



M + ^^k)Srn-nS + (4' " )'5r, 



(29) 



The Hamiltonian part of Eq. (|29p is explicitly written down in a tridiagonal matrix form as 















Cke^''' 

















Cke-'"" 





CkC''^" 












(30) 



r 



As remarked in Sec. [iTl each component of the Floquet 
matrix [Eq. ([50)1 ] represents a probability amplitude of 
a transition from one Floquet mode to another. Note 
that for the case of the dc field, Hamiltonian ([50)1 has 
no diagonal components, which means that the electrons 



cannot stay stationary but are always excited by the field. 
We also note that the off diagonal components do not 
depend on Q which is proportional to the strength of the 
field. The dependence of fl is only taken into account 
through the quasienergy part of Eq. ([29]) . 
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FIG. 1: The coefficients ae of the delta functions in the local 
spectral function of the noninteracting electrons on the hy- 
percubic lattice with the dc field fl = 0.25 are plotted by the 
circles. The delta functions with the spacing Q are schemati- 
cally shown by the solid lines. The broken line is a guide for 
the eyes. 



2. Hypercubic lattice in an ac field 

Let us move on to the case of the ac field on the hy- 
percubic lattice. The vector potential is defined by 

eA(t)a = AsinUt, (31) 

where is the frequency of the ac field, and 

eEa 



A = -- 



(32) 



is its amplitude divided by the frequency. Although we 
use the symbol 17, this should not be confused with Vl (the 



Bloch frequency) introduced in Eq. ([24]) for the dc field. 
Following Eq. we derive the Floquet representation 
of Green's function as 



(G«V„(c.)=5]. 



1 



dx r dy 



— TT <J —TT 



2n 2tt 



X exp 



dz {e^''[cos(Asinz) - Jq{A)] 



sin(Asinz)} 



(33) 



Taking the imaginary part of Eq. (|33p and integrating 
over k with the JDOS [Eq. (|E2|) ]. we obtain the local 
spectral function. We depict it for several A and = 1 
in Fig. [21 One can see that the spectral function has nar- 
row peaks at uj = nfl {n = 0, ±1, ±2, . . . ) known as the 
dynamical Wannier- Stark ladder. The peaks at u = ±f2 
correspond to one-photon absorption or emission, and 
the peaks at lu — ±2f2 correspond to two-photon one, 
etc. The width of each peak, or the effective hopping, is 
renormalized by the zeroth-order Bessel function Jo{A) 
as seen in Eq. ([33]) . and even vanishes making the elec- 
trons completely localized when Jo(^) = 0. In Fig. [2] 
one finds that the widths of the peaks shrink as A ap- 
proaches the first zero (z = 2.40483 ••• ) of Jo{A) until 
finally the peaks become the delta functions. This scal- 
ing has been known as dynamical localization since the 
proposal by Dunlap and Kcnkre.^'^ 

The inverse of Green's function in the ac field can be 
calculated via Eq. (flSl) . The result is 



J 



(G^° )„,n{oj) ^ {co + nn + n + iri)5^ 



e-kJ-m-niA) m 
le'^Jm-niA) m 



n: even 
n: odd 



The explicit Floquet matrix form of the Hamiltonian part of Eq. ([M]) reads 



efeJo(A) 


iefeJi(A) 


ekJ2{A) 


iekMA) 


€kM^) 


-ilkJi{A) 




ilkJiiA) 


ekMA) 


iekMA) 


f-kJ2{A) 


-iekJi{A) 


ekMA) 


iekMA) 


ekMA) 


-iikJ^iA) 


ekJ2{A) 


-iefcJi(A) 


ekMA) 


iikMA) 


ekMA) 


-iekJsiA) 


ekMA) 


-iekMA) 


ekMA) 
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(34) 



(35) 



where we omit the label "sc" attached to ek and efc for 
simplicity. Taking the dc limit — > requires a special 



care because of definition ([32)1 which is singular at Q — 
0. Therefore a quantity calculated for a system in the 
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4ii 





1 A = 






1 





















fee 



Q-1 



= ^ ^ COS fca COS A:/3 . (36) 



2 ^2/5=1 

In the infinite-dimensional limit {d oo) the dispersion 
of the fee lattiee e^'^ 
[Eq. dUl)] through^ 

1 



of the fee lattiee e^'^ is related to that for the sc lattice 



(37) 



FIG. 2: (Color online) The local spectral functions of the 
noninteracting electrons on the hypercubic lattice coupled to 
the ac field with = 1 and A = 0, 1.5, 2, 2.2, and 2.3. 



Using Eqs. and ([57]) . we derive the inverse of Green's 
function for the infinite-dimensional fee lattice in the de 
field, 



presence of the ac field with a finite frequency ^ 
does not necessarily reproduce the result calculated for a 
system with the dc field discussed in Sec. IIV C II 

Let us examine the physical meaning of the Floquet 
matrix [Eq. (|35p ]. The (rn,n) component of Hamilto- 
nian ([35|) is proportional to Jm-n{A). Since Jm~n{A) cx 
j[\m-n\ j-j- ^ sufficiently small, the transition proba- 
bility m ^ n is proportional to E'^^"^~^K For m < n, 
the process corresponds to stimulated absorption, while 
it corresponds to stimulated emission for m > n. The 
process of spontaneous emission is not included since we 
assume that the electric field is classical so that there 
is no quantum fiuctuation of photon numbers. This as- 
sumption is appropriate as long as the intensity of the 
electric field considered here be so strong as a pulsed 
laser. 



D. Application to other lattices 

So far we have assumed that the vector potential A 
points to the specific direction (1,1,...,!) in the hyper- 
cubic lattice. We can more generally calculate the inverse 
of the noninteracting retarded Green's function G^" 
by making use of formula (jlSp for arbitrary lattiee struc- 
tures and the vector potentials. The disadvantage of such 
a general ease is that, since the fc dependence of 
is not so simple as to be only through e^^ and ej^, the 
integral over fc becomes computationally heavier. 

Here we note that there is a group of lattice models 
that give a simple fc-dependenee of G^''^ in infinite 
dimensions. Among them are the face-centered cubic 
(fee) (Sec. IIV D II) and body-centered cubic (bee) 
(Sec. IIV D 21) lattices with A parallel to (1, 1, . . . , 1) in 
infinite dimensions. It is instructive to give examples 
other than the hypercubic lattice. 

1. fee lattice in infinite dimensions 

The energy dispersion of the fee lattice generalized to 
arbitrary d dimensions (d > 2) is given by 



(G™ )nin{^) = icj + nQ + fl + iri)S„,n 



+ 2{Cl - 1)(5,„„] 



(38) 



We see that the Hamiltonian part of the inverse of 
Green's function [the second term on the rhs of Eq. (|38p ] 
is written in a pentadiagonal matrix form. In the same 
way we can obtain Green's function on the ac field via 
Eq. ([HI). The resuh reads 



(G™ )mn{^) = (o; + nfi + /i + iri)S„in 

1 fClcOsi2ek)Jm-n{'2A) + {Cl-l)S,nn m 



2 HCfcSin(26'fc)J„_„(2A) 



m — n: odd 
(39) 



We notice that the factors Jm_„(2^) and Jm-n(O) = Smn 
appear in Eq. ((39|) [while Jm-niA) appears on the sc lat- 
tice; see Eq. ([M)) ]. Equations ([55)) and ([55]) indicate that 
Green's functions depend on k only via the two functions 
Efc [Eq. ([HI)] and e"^ [Eq. so that we can integrate 

over fc using the JDOS [Eq. ([52l) ]. 



2. bee lattiee in infinite dimensions 

The bee lattice in d dimensions {d > 3) is defined by 
the dispersion relation. 



bcc 



g d Q-1/3-1 

. > > > COS k„ cos kfj cos 

2./d{d^T){d^) hhh 

We can take the limit d — > oo in the same way as in the 
case of the fee lattice, and the dispersion converges to 



(41) 



from which we can derive Green's function in the dc field. 
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(Gf ^),nnH ^{u; + nn + ^i + Z7?)<5„„ - — e'("-")'''=Cfe[Cfc('5m-«,3 + <5m-„,-3) + 3(Cl - 2)(<5„_„,i + <5„._„,_i)]. 

(42) 

In this case the Hamihonian part of the hiverse of Green's function [the second term on the rhs of Eq. p2|) ] becomes 
a heptadiagonal matrix. Similarly, Green's function in the ac field is written as 

rnRo^h ( \_( , oi , • \a 1 /C/=[Cfc cos(36lfc) J,„_„(3A) + 3(Cfe - 2) cosSfc J,„_„(A)] m - n: even 

6 [iCfcKfcSm(36'fc)J„i_„(3yi) + 3(Cfc - 2)sm6'fcJ,„_„(A)J m - n: odd 

(43) 



where the factor J„i_„(374) newly appears besides 
Jm-n(^)- Again, Green's functions depend on k only 
through e^j^ and e^^^ , which makes an integral over k com- 
putationally quite efficient with the JDOS [Eq. 



DIAGONALIZATION OF FLOQUET 
MATRICES 



Here we examine how to diagonalize a Floquet matrix 
{lo -^n^)5mn — Hmn appearing in the original Schrodinger 
Eq. ([U for the noninteracting electrons. We have shown 
in Sec. IIV Al that the Floquet matrix form of Green's 
function is diagonalized into Qk{^) by the unitary trans- 
formation Afc [see Eq. ((T7)) ]. In Sec. lIVBl we have men- 
tioned that the inverse of Green's function is equivalent 
to the Floquet matrix (w -I- nfl + + iri)Smn — Hmn ■ Com- 
bining these two facts, we identify the eigenvalues and the 
eigenvectors of the Floquet matrix. We summarize the 
statement below. 

The eigenvalues of a Floquet matrix H„in ~ nflSmn for 
a single-band noninteracting system subject to a homoge- 
neous electric field periodic in time are 



nn (n = 0,±l,±2,...), (44) 



and for each n, the corresponding eigenvector is given by 



"" = (Afc)„„ (m = 0,±l,±2,. 



(45) 



This completely solves Floquet matrix problems for a 
single-band system of noninteracting electrons. We note 
that the fact that Eq. (|^^ gives the eigenvalues is es- 
sentially a consequence of Shirley's relation [Eq. (5) of 
Ref. [lal- Since (efc)oT is a dynamical phase, the above 
statement indicates the absence of an additional geomet- 
rical phase (see Appendix [P]) . 

From Eqs. ([TT]) and the noninteracting Green's 

function can be written with the Floquet states as 



to + £Q -\- fj, — (efc)o + it] 



(46) 



One should note that the theorem cannot be applied to 
a multiband system, where contributions of a geometri- 
cal phase may survive and interband transition could be 
caused by the electric field. 

The theorem provides a unified description of period- 
ically driven systems: the original band structure efc is 
renormalized into the time-averaged one (efc)o due to the 
field, and the renormalized band splits into its replicas 
with the spacing fl. On the (hyper)cubic lattice, the 
dc field changes the band into the Wannier-Stark lad- 
der with the infinitesimal bandwidth [Fig. ^eC)] , while in 
the ac field, the bandwidth scales with the factor Jo{A) 
[Fig.Hb)]. 

To see how the theorem actually works, let us apply 
it to the d-dimensional fee lattice model in the dc field 
as an example. For convenience, we restrict our discus- 
sion in the limit d oo. The Floquet matrix form of 
the Hamiltonian is given by Eq. (j38p . According to the 
theorem, its eigenvalues are equal to the diagonal com- 
ponents: uj -\- nfl — ^{Ck ~ !)■ This means that the dc 
field modifies the energy dispersion from e|"^ [Eq. ((37)) ] to 
e^^^ = ^(Cfc"!)- The latter dispersion is equivalent to the 
one for a (d — l)-dimensional hyperplane perpendicular 
to(l,l,...,l), the direction of the electric field. We can 
interpret this fact as follows: the dc field makes the elec- 
trons localize along the direction of the field due to the 
Bloch oscillation, but the electrons are free to move along 
the directions perpendicular to the field because along 
those directions the force of the dc field does not act on 
the electrons. As a result, the motion of the electrons 
is confined on the hyperplane, and the energy dispersion 
becomes e^^'^. 

Another example is provided by an ac field. If one sees 
the Floquet matrix form of the Hamiltonian for the sys- 
tem in an ac field on the fee [Eq. or bcc [Eq. (^5)) ] 
lattices, one finds that various kinds of band renormaliza- 
tion besides the Jq (^) scaling on the simple cubic lattice 
are implied by the theorem. 

Finally, let us consider electrons in a parabolic band 
(e.g., a conduction-band bottom in a semiconductor) 
with efc = k^/2m*. If we apply an ac field on it, the 
dispersion is renormalized into (efc)o = ek -\- e'^ E'^ / Am* fl^ 
[Fig. [3Ic)]- This results from the dynamical Franz- 
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FIG. 3: The band renormalization €k (efc)o on the (hy- 
per) cubic lattice in (a) the dc field and (b) the ac field and 
(c) on a parabolic band in the ac field. The solid line repre- 
sents Efe, while the dashed line respresents (efe)o. The circles 
are occupied states moving along the arrows. 



Keldysh effect,— which is, in the present formaUsm, 
naturally understood through the diagonalization picture 
of the Floquet matrix. 



VI. DYNAMICAL MEAN-FIELD THEORY 
WITH THE FLOQUET-GREEN FUNCTION 
METHOD 

Now we are in position to combine FGFM with DMFT 
for treating interacting systems in external fields. At 
the basis of DMFT lies the fact that a lattice prob- 
lem of correlated many-body systems can be approx- 
imately mapped to a problem of an impurity embed- 
ded into the environment of an effective medium when 
one ignores spatial fluctuations but takes fully into ac- 
count on-site dynamical correlationfi2i21 The mapping 
is given as follows: let Z = f[dci][dcl] e^'^I'^''^^! be 
a partition function in terms of the original action, 

Jdtjdt' 4 it)Gl'\t, t')c, (f) + ^int h , 4] : 

where the interaction term is assumed to be a sum 
of the local terms. Integrating out each site's de- 



grees of freedom except for a representative site i = 
o, we have the local partition function, Ziad'^o] = 
J[dCo][dcl] e*'^i°=[^°''^<'l. Here the local action reads ^loc — 
Jdtjdt'cl{t)%~\t,t')c,it') + S,r.t[co,cl] [%{t,t'): the 
Weiss function]. If one ignores the nonlocal fluctua- 
tions, Z and Zioc give a common site-diagonal self-energy, 
T,ij{t, t') — 5ijYj{t, t'). This fact enables us to build a set 
of self-consistent closed equations, which can be solved 
with an iterative numerical calculation. Although ne- 
glecting the spatial fluctuation is generally an approx- 
imation, the nonlocal corrections become rigorously ir- 
relevant in the limit of infinite dimensions, where the 
hopping parameter is scaled ast = t* /2y/d (t*: fixed)i^ 
DMFT is also applicable to nonequilibrium systems as 
recently studied^^i^i^ based on nonequilibrium Green's- 
function formalism. Similar to the equilibrium case, one 
has self-consistent equations for Green's function and the 
self-energy. Now, the present proposal is that if the 
driven system is periodic in time, one can rewrite the 
equations in the Floquet matrix form. 



(Gloc)mn('^) 

(Gloc)ni"(^) 
(Gfc^)nin(^) 
(Gloc)mn('^) 



5Zi,,[%] 

(G^ )mn{u!) — Y,mn{tij), 

^(Gfc)m„(a;). 



(47) 

(48) 
(49) 
(50) 



To solve Eqs. (|T7)) - (15(71) self-consistently, we first in- 
put the inverse of the noninteracting Green's function 
given by Eq. HH) into Eq. ([49]). After the initial self- 
energy is properly chosen, the calculation is iterated un- 
til Green's function converges. For the lattices discussed 
in Sees. IIV Cl and HV Dl the integral over k in Eq. ([501) is 
performed via the JDOS [Eq. ((22)) ]. 

As remarked in Sec. IIIIl the size of the Floquet matrix 
that needs to be taken in a calculation is usually small 
(~ 5-30, depending on $1), which, with the analytic ex- 
pression of the inverse of Green's function (jl8[) . makes 
our computational costs dramatically small. 



VII. GAUGE-INVARIANT GREEN'S 
FUNCTION 



Before applying our method to a model, we exam- 
ine the gauge invariance of Green's function. Let us 
write the coordinates x"^ = {t, r) and the vector po- 
tential A"^ — (0, A) in the four-vector form. The 
gauge transformation, Ai,{x) A^{x) + d^xi.'^)-, Puts 
a phase factor to the creation and the annihilation oper- 
ators as c^{x) e~^'^^''^^ {x) and c{x) e*'^'^^^-'c(a;). 
Accordingly Green's function changes as G{x,x') 
g«e[x(2:)-x(2: 3;')^ i.e., the usual Green's function 

is not gauge invariant. It is known^i^^ that one can 
make Green's function gauge invariant with an additional 
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phase factor as 

G{x, x') = exp ( -^ / dy" eA,{y) ] G{x, x'). (51) 



G depends on the path of the Hnc integral m the expo- 
nent. Here we adopt the conventional straight line con- 
necting X with x' as the path of the integral. Suppose 
that Green's function depends on k only through and 
ej^. Then in the temporal gauge (0 = 0) the Wigner rep- 
resentation of the modified Green's function G becomes 



dt 



rcl 



X c 

y.Gn{^,0 + j ^ d\ eA(tav + A^rd), t^'^ : 

(52) 



where C, and 9 are defined in Eqs. ([26|1 and (p7|) . Note 
that G is calculated by shifting the variable 9 in the orig- 
inal Green's function. This suggests that Green's func- 
tion integrated in terms of 9 is definitely gauge invari- 
ant, so that the local Green's function G'fc(w) is also 
gauge invariant."23 In the following, we derive the gauge- 
invariant Green's function G for the dc field in Sec. lVIlXl 
and the ac field in Sec. IVII Bl 



A. dc field 

To obtain G for the dc field, we first note that the 
Hamiltonian in the dc field [Eq. (|23p ] has the time trans- 
lation symmetry. If one makes a time translation t — > 
t + 5t, the vector potential changes as A{t) A{t) + ^lSt. 
Since the change can be absorbed by the gauge trans- 
formation with X = ~^^^'^'i=i^^i the Hamiltonian is 
invariant against time translation. Then we assume that 
in the long-time limit after the dc field is switched on the 
retarded Green's function becomes independent of the 
initial correlations. This assumption seems to be valid^ 
as numerically checked in Sec. IVIII Al As a result, a 
gauge-invariant quantity that is calculated from the re- 
tarded Green's function should be independent of the av- 
erage time tav For instance, the local Green's function 
in the Wigner representation (Gj^j.)„(w), which is gauge 
invariant as shown above, vanishes for n ^ 0, so that 
Gj^j.(i,t') does not depend on tav In the same way the 
self-energy S^(w) also vanishes for n ^ 0. 

Since the Floquet-represented self-energy (w) is di- 
agonal due to the symmetry, we can identify the 9 de- 
pendence of the Floquet representation of the retarded 
Green's function as 



G«„(C,0,c.)=e'("-")^G^„(C,0 = O,c.). 



Using Eq. ([53]) . one can evaluate the gauge- invariant 
Green's function (|5^ as 



G^(C,0,c^) =5™,o^G^'(C,e,c^), (54) 

n 

where we can see that every mode of Green's function 
equally contributes to G^. 



B. ac field 

For the ac field, which is one of the key questions in the 
present paper. Green's function has a more complicated 
9 dependence. To evaluate the gauge-invariant Green's 
function (|52p . here we expand the original Green's func- 
tion with respect to eA{t) in a Taylor series: G^{(^,9 + 
JdX eA,u;) = EiUJdX eAf a^G«(C, 0, o.). Then 
the Wigner representation of the gauge-invariant Green's 
function is expressed as 



in ■ •' 

xxi^L„r(^)(^(c.-c.')), (55) 



where 



(2 



m — n.2r—i-i 



(56) 



and 



' , 27r \ X 



3T)! ^ C) + ^ 2r)^-i0(fc + £-2r) 

(57) 



2^(^-1) 



r=0 



for i > 1. When £ = 0, we have y(^)(fc) = S{k). The 
Floquet representation of Eq. ((55)) reads 



2 /A 



G^^{C9,Lu)^GtACd,^)+J2 



1=1 



r=0 



f.n/2 

+ 

A;=0 "'-^V2 fe^iy 

(58) 



(53) 



The derivative with respect to 9 in Eqs. (|55p and ([55]) 
is simplified when the system is on the hypercubic lat- 
tice. The Floquet representation of the noninteracting 
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Hamiltonian H is then given by Eq. ((35)) . and the £th 
derivative of the Floquet representation of the retarded 
Green's function can be calculated for every I via the 
recurrence relations 

deH = iJ, (59) 
deH = -H, (60) 
dgG'^ = -G^{-dgH)G^ = G^HG^. (61) 

Employing Eqs. ((59l) - ((6T|l with relation ((58)l . one can nu- 
merically evaluate the gauge-invariant retarded Green's 
function G^ for the ac field. 



VIII. APPLICATION TO THE 
FALICOV-KIMBALL MODEL 

To test the ability of the present method for treat- 
ing many-body systems, we apply it to the spinless FK 
model, for which the Hamiltonian is 

H = -Y. t^A^o + uY^cUflh. (62) 

ij i 

Here fi (//) annihilates (creates) a locahzed electron, 
and [/ is a coupling constant. It is known that the FK 
model exhibits a metal-insulator transition in infinite di- 
mensions from DMFT calculations, where the possibility 
of charge-density wave phases is ignored.— The critical 
value of U for the transition is known to be on the 
hypercubic lattice at half filling. The insulating phase 
is Mott-like, which means that the insulating state orig- 
inates from the electron correlation. 

What characterizes the FK model is that it has an 
exact solution for the impurity problem [Eq. (|47p ] within 
DMFT, even out of equilibriumj^ The solution for the 
retarded Green's function is 

G^iLu) = wo %^(c.) + w, [%'''\co) - U]-\ (63) 

where wi is the filling of the / electrons, and wq — 1 — wi. 
Note that the retarded component of Green's function de- 
couples to the Keldysh component. Here we concentrate 
on the retarded Green's function, and calculate the lo- 
cal spectral function A„(a;) = —^lm{G^^)nii^) (which is 
gauge invariant as explained in Sec. IVII[1 and the gauge- 
invariant spectral function An{k,ijj) — —^lm{G^')n{^^) 
under the assumption that the initial correlations are ir- 
relevant to the retarded Green's function. 



A. Falicov-Kimball model in the dc field 

We first present the results for the dc field. We start 
with noting that the integral over 9 can be performed 
analytically due to relation (|53p . and that Gj^^., ^^q^, and 
in the Floquet representation are all diagonal as men- 
tioned in Sec. I VII XI which simplifies our calculation. All 



U=0.5 U=2.0 




-2-10 1 2 -3-2-10123 



FIG. 4; The local spectral function Ao{uj) for the FK model 
coupled to the dc field on the hypercubic lattice at half filling. 



the matrices we have to invert numerically are tridiagonal 
owing to Eq. (^5)) . 

In Fig. [4] we illustrate the local spectral function Ao{uj) 
for various values of U and f2 on the hypercubic lattice. 
The size of the Floquet matrices that we choose is typi- 
cally 9-13. Convergence is achieved after typically 10-30 
iterations, where the calculation is quite stable over the 
parameter space considered here. We can see that the 
present result in Figs. |4Ub)-md) obtained in the Floquet 
method agrees well with the previous results, where 
the nonequilibrium DMFT is employed. This suggests 
that our assumption of the irrelevance of the initial cor- 
relations to the retarded Green's function is valid. In 
our results, the spectral function is positive definite, and 
satisfies the sum rule for the zeroth spectral moment^^ 
as in equilibrium. Therefore we can safely interpret the 
quantity A{ui) as the spectrum of the system even out of 
equilibrium. 

More interesting case is Figs. |3Ue)-|3l^h), where we can 
observe how a Mott-like insulator {U — 2) is driven into 
a metallic state by the dc field. Namely, while there is 
a clear band gap between the upper and lower bands in 
equilibrium [Fig. I^e) , = 0] , the gap disappears as the 
dc field is increased, where the spectral weight around 
uj = develops. Hence our calculation captures the Mott- 
like insulator-to-metal transition induced by a static elec- 
tric field. In the strong dc field region [Figs. [3Kb)-|3Kd), 
Sfg), and[Hh)], we find complicated structures with the 
spacing SI. We can attribute these to the Wannier-Stark 
ladder (mentioned in Sec. IIV C 1[) . which grows with the 
field intensity E (x ft [see Eq. The Wannier-Stark 

structure interferes with the original spectrum that com- 
prises two bands with the spacing ?7, producing a char- 
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FIG. 5: The zeroth mode of the local spectral function Aq{ijj) 
of the FK model coupled to the ac field {Q, = 1) on the hy- 
percubic lattice at half filling for U — 0.5. 



acteristic interference pattern. 



B. Falicov-Kimball model in the ac field 

We move on to the ac field. While the spectral function 
is time independent for the dc field, i.e., the non-zeroth 
modes of An{u!) vanish, this is no longer the case for the 
ac field. Since our interest resides in the time-averaged 
spectral function J dti^v ^(w,tav) = Aq{uj), we concen- 
trate on the zeroth mode of the spectral function. Unlike 
the dc case the integral over 9 is nontrivial, which has 
to be calculated numerically. In Figs. [5][51 we depict the 
local spectral function Ao{uj) on the hypercubic lattice 
at half filling with the frequency of the ac field f2 = 1. 
The efficiency of the convergence and the stability of the 
calculation are similar to the dc case. 

In the metallic region (Fig.[5]for U = 0.5 and Fig.[6]for 
U = 1.3), one can see how the metallic spectrum of the 
system is deformed by the ac field. Namely, the width 
of the band shrinks with the intensity of the field. It 
can even goes to zero when A coincides with a zero of 
Jo{A), which makes interacting electrons localize. This 
is quite similar to the noninteracting case as examined 
in Sec. IIV C 21 Hence we have the dynamical localization 
in interacting electron systems. Note that the scaling of 
the band width with Jo{A) is a nonlinear effect of the 

ac field, as evident from Jo{A) = 1 - {A/2)^ H . The 

difference between noninteracting and interacting cases 
is that each peak in the dynamical Wannier-Stark ladder 
at w = nil {n = 0, ±1, ±2, . . . ) splits into two with the 
spacing U due to the correlation effect. This can clearly 
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FIG. 6: The zeroth mode of the local spectral function Ao{uj) 
of the FK model coupled to the ac field {fl — 1) on the hy- 
percubic lattice at half filling for U = 1.3. 



be seen in Fig. [5^c). 

In the insulating region (Fig. [7] for U — 2.2 and Fig. [5] 
for U = 3.8), on the other hand, we do observe the 
ac-field driven transition from the Mott-like insulating 
state to a metallic state. In equilibrium [Figs. [7^a) and 
[SJ^a)] the system has a gap between the upper and lower 
bands. When the ac field is switched on, the gap collapses 
[Fig. Efb) and[5](b)], and a spectral weight grows in the 
midgap region around uj — 0. If we compare Figs. Wljo) 
and[5](b), we can see that the larger the band gap, the 
smaller the midgap weight. From these results, we see 
that metallic states appear in the insulating system of 
correlated electrons in the intense ac field. As the inten- 
sity of the ac field is further increased, the system plunges 
into the dynamical localization regime, where the band 
width starts to scale with Jo{A). 

To characterize the metallic state, we have calculated 
the momentum resolved spectral function Ao^kjUj). For 
clarity we take the simple cubic lattice with the JDOS 
[Eq. (|E5[) ] in three dimensions. As a key result in the 
present paper, we plot the zeroth mode of the spectral 
function Ao{k,oj) and Ao{llj) for the frequency 17 = 1 at 
half filling in Figs. and [TUl where we take k = fc(l, 1, 1). 
One can check that the result respects the particle-hole 
symmetry. As discussed in Sec. IVIIBl numerical eval- 
uation of Ao(fc,a;) is done in a perturbative way. Thus 
we can obtain reliable results only in a weak intensity re- 
gion. Although we use a perturbation in A, higher-order 
contributions are included in our calculations. To obtain 
the results in Figs. [5] and 1101 the summation over £ in 
Eq. ([55]) is performed up to ^ = 5. We have checked that 
the expansion in terms of A converges for A < 0.6. A 
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FIG. 7: The zeroth mode of the local spectral function Ao{u}) 
of the FK model coupled to the ac field {Q — 1) on the hy- 
percubic lattice at half filling for U — 2.2. 



FIG. 8: The zeroth mode of the local spectral function Ao{uj) 
of the FK model coupled to the ac field (57 = 1) on the hy- 
percubic lattice at half filling for U = 3.8. 



calculation tends to be unstable when U is small (< 1) 
where the system is in a metallic state. When we an- 
alyze the results, we have to be careful with the sign 
of the spectral function Ao{k,uj). Although the local 
spectral function Ao{lu) is positive definite, Ao{k,Lu) is 
not so in general. While we notice there are some re- 
gions where Ao(fe,a;) becomes negative, the quantity is 
mostly positive for A < 0.6 and U >2. The result should 
be supported by other gauge-invariant quantities such as 
the current or the optical conductivity, which is a future 
problem. 

First, let us see Figs. |9l[a)-[9l[c). These are the spectra 
of the system in equilibrium (A = 0). One can see how 
the metallic bands [Fig. O^a)] change into the insulating 
ones [Figs. Elb) and [He)] with a finite gap appearing 
with U. In the insulating state, the upper band is almost 
a replica of the lower one shifted upward by U, which 
is characteristic of the FK model. When the ac field is 
turned on in Figs. fTUfa^-fTUTc) {A — 0.6), we can see how 
the ac field generates a new photoinduced band structure. 
In Fig. fTOl b). we can observe that a new band appears in 
the midgap region. This band is created by the electrons 
that absorb or emit one photon with the energy 17 = 1, 
that is, the photoinduced band is a replica of the original 
lower and upper bands shifted by 51. If we assume that 
the states are occupied up to w = as in equilibrium, 
the electrons in the induced band around a; = play 
a role of carriers, making the system metallic. When 
the interaction U is strong enough [Fig. [TOFc). U = 3.8], 
the metallic band does not appear. Instead, side bands 
appear near the original bands with the spacing in the 
midgap region. Again, the electrons in the side bands 



consist of the electrons absorbing or emitting one photon 
with the energy fl. Since il is much smaller than U here, 
the electrons cannot reach the region around lo = with 
a one-photon process. Therefore the system remains to 
be insulating with the finite gap in the ac field. As well 
as the case of the dc field, the side band pattern with 
the spacing Q interferes with the original band structure 
with the spacing U, yielding complicated band patterns. 



C. Relevance to experiments 

Finally we mention the relevance of the present re- 
sults to experiments. For the dc field, the intensity re- 
quired for the effect considered here is ^ 10^"^"^ V/m for 
a ^ 10°"^A, which is too strong to be realistic. However, 
in the case of the ac field, the required intensity is A ~ 1 
[see Eq. ([5^ ] for the dimensionless quantity, which trans- 
lates to £; ~ 10^-1° V/m (i.e., the intensity of ~ 10""" 
W/cm^) for ~ 1 eV (visible fight). For a smaller fl the 
required intensity becomes smaller. Since the intensity 
of pulsed laser available with recent advances in optical 
techniques reaches such magnitudeSf^I it should be pos- 
sible to observe the nonlinear effects predicted in this 
paper in experiments. One problem is that when the 
intensity goes beyond ~ 10"'^^ W/cm^, atoms begin to 
be ionized and evaporated. To make the required field 
intensity smaller, we can take systems with large lattice 
constants, such as the zeolites loaded with guest atoms.— 
As an entirely different class of systems, we can consider 
cold atoms in optical lattices^^S which may be an interest- 
ing playing ground for the effects examined in the present 
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FIG. 9: (Color online) The zeroth mode of the gauge-invariant 
spectral function Ao{k,u!) with k — fc(l, 1, 1) (the density 
plots) and the local spectral function A(){uj) (the line plots) 
of the FK model coupled to the ac field (fi = 1) on the simple 
cubic lattice at half filling for A = in units of t* . The color 
bars on the right side represent the correspondence between 
the colors and the values of Ao{k,Lj). 




FIG. 10: (Color online) The zeroth mode of the gauge- 
invariant spectral function Ao{k,ui) with k = A:(l, 1, 1) (the 
density plots) and the local spectral function Ao{uj) (the line 
plots) of the FK model coupled to the ac field {fl — 1) on the 
simple cubic lattice at half filling for A = 0.6 in units of t* . 
The color bars on the right side represent the correspondence 
between the colors and the values of Ao{k,u!). 



paper. 



IX. CONCLUSION 

We have developed a theoretical method to formulate 
photoinduced phenomena in correlated electron systems. 
The method incorporates FGFM into DMFT, which can 
fully take into account both the electron correlation effect 
and the nonlinear electric- field effect. We have applied 
the method to the Falicov-Kimball model in ac fields to 
calculate the gauge-invariant spectral functions. We find 



peculiar photoinduced band structures, which arise from 
the nonlinear effect of the electric field. In particular, we 
find a metallic state in the midgap region of the Mott-like 
insulator induced by the ac field. In the calculation we 
have utilized a theorem, found here, that identifies eigen- 
values and eigenvectors of a Floquet matrix for single- 
band noninteracting electrons. 

We believe that our approach has a potential ability 
to treat, not only the FK model considered in the pa- 
per, but also a wide class of models such as the Hubbard 
model. There are some future problems: one is to calcu- 
late the Keldysh component of Green's function . We 
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need the Keldysh Green's function to compute, e.g., the 
current or the optical conductivity which has information 
of the transport properties of the system. An apphcation 
to the Hubbard model is also desirable. Experimentally, 
the relaxation of photoinduced states after the ac field is 
switched off is also an important phenomenon. This is 
theoretically interesting as well, for which further devel- 
opments on the nonequilibrium DMFT would be neces- 
sary. 
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APPENDIX A: MULTIPLICATION RULE FOR 
THE FLOQUET MATRICES 

Here we show that Floquet matrices obey the multi- 
plication rule in the linear algebra. Let us prepare two 
functions A{t, t') and B{t, t') which satisfy the periodicity 
condition: A{t + T,t' + t) = A{t, t') (and so does B). We 
write the following integral in the Wigner representation: 



E 



2tt 

dhJ 

2^ 



-i(w-i'n/2){t-t')-i(i+t)n(t+t')/2^^l^^-^Q^^ 



£ + i' . 



-n 



r 



Thus we have 



{AB)k{Lu)^ Ai(uj+^-n^Bi,[uj~^n 



Let us take an integer k' satisfying the condition k' 
k (mod 2). Replacing oj with oj + k'il/2 gives 



E Ae[uj + ^^n 

/ k' — I 
X Bfluj H — -fl 



If we write Eq. (jAip in the Floquet representation fol- 
lowing its definition ([7]), we arrive at the conclusion, 



(^■B)(fe+fe')/2,(fe'-fe)/2(^) 
= E '^{k+k')/2Xk+k')/2-e{^)B(^k+k')/2-e,{k'-k)/2{'^) 

e 

or a more transparent expression, 

mm' ij^^BrfYi'ni}^^-; 



to k = k' (mod 2). This ensures that we can apply 
the usual multiplication rule of a matrix to every 
Floquet-represented function. 



APPENDIX B: DERIVATION OF THE FLOQUET 
REPRESENTATION OF THE 
NONINTERACTING GREEN'S FUNCTION 



Here we derive the Floquet representation [Eq. (|13p ] 

(Al) 

of Green's function. Let us start with Eq. We 
find that the argument of the exponential in Eq. (fT2l) is 
not invariant under discrete translation against tj-d. To 
somehow make it invariant under such a translation, we 
rewrite Eq. p2)) into 



/ /.tav+t„l/2 ^ 

X exp I -i I dt" [€k-eA(t") - (efe)o] 



tav-irel/2 



(Bl) 



™ Here (efc)o is defined in Eq. (|14p. Now that the argu- 

where m — {k + k')/2, n — {k' — k)/2, and ment of the exponential in Eq. (|B1[) is periodic in t^ci 
m' = {k -\- k')/2 — all of which are integers due with the period 2t and in tav with the period r, we can 
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insert the factors «^f^*roi/2_L dt'^^y Q^^^'^trd/^ and where 77 is an infinitesimal positive constant, we can per- 

^-^mat^^ i ^^^^ ^^mat'S Eq. jBll) . Then, Wigner transformation of Eq. m as 

with the Fourier transformed expression of the step func- 
tion, 

0{tr-ci) cLu' — (B2) 

ZTTi J LlJ + iri 

I 



^7 2TIUJ' + If^j tJ_^/2 



2t 



dt' 



dC e^^^^-Z^+^^f^^v exp dt" [ek-eA(t") ~ (efc)o] ) • (B3) 



-T/2 



*av-*;cl/2 



In order to make our notation clearer, we change the 
integral variables as — x', fit^v = y' ■ After some 

calculations, we obtain 



1 



^ w - en/2 + ^ - (efc)o + ir] 



mod 2 



X / /" iix'+iny' 



X exp 



y'+x' 
^ Jy'-x' 



dz [ek-eA(z/n) — (efc)o]^ 
(B4) 



In the above we have used the fact that the inte- 
gral in Eq. (|B4p equals zero when £ ^ n (mod 2) 
since JJdx'dy' = + + (//„ + //iv) = [1 + 

(— l)^+"](JJj + JJjj), where the Roman numerals repre- 
sent the ranges of the integral defined in Fig.[TT] We fur- 
ther change the integral variables: x' + y' ^ x, x' — y' = 
—y. Here we have to be careful with the range of the 
integral. As shown in Fig. [Til we change the range of 
the integral from J JJdx'dy' to ^ J JJdx'dy', which is 

equal to i /^2ir dx /^27r ^y timcs i coming from the Ja- 
cobian. Then we change the range of the integral again: 
i X i /^27r d^ /-27r ^y — /T^ dx J^^ dy . Notc that the two 
i factors are canceled out by the change of the range of 
the integral. As a consequence, we arrive at the general 
Wigner representation of Green's function. 



(Gj?")n(^) 



E 



^ uj ~ m/2 + fi — {ek)o + if] 



mod 2 



dx 
2^ 



dy 
2^' 



,i{i+n)x /2-i{e-n)y /2 



X exp 



dz [ek_ 



eA{z/n) 



(efc)o] 



(B5) 
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FIG. 11: The range of the integral: each Roman numeral 
denotes the corresponding shaded region, the symbols O and 
□ used in the text denote the regions surrounded by the solid 
line and the broken line, respectively. 



Next, let us move on to the Floquet representation. 
Following the definition of the Floquet representation 
[Eq. (g])], we have 



E 



1 



mod 2 



uj + (m + n — ^)ri/2 -I- ^ — (efc)o + if] 
dx r dy 



27r 



J.{l+m-n)x l2-i{t-m+n)y j2 



2-K 



X exp 



dz [efe_eA(z/n) - (efc)o_ 



In the above we can replace m ^ n — i with 2i due to 
m-fn — ^ = TO — n — £=0 (mod 2) , which gives Eq. ([Te 



17 



APPENDIX C: UNITARITY OF Afc 



We prove that defined by Eq. (|T5)) is a unitary 
matrix for any e^. and A{t) as 



BA.).-..(Ai)'..^/'J/'Ji: 



^i{rfix~ny)~i£[x~y) 



X exp J dz [ek-eAiz/n) - (efc)o] 

— / dy S{x - y) e'('"^-«f) 

exp (^^-^ J [efc-eA(^/n) - (efc)o] 



dx 



— O77 



Since is nothing but a set of Floquet state vectors 
[Eq. (j45|) ]. this unitarity relation indicates that the Flo- 
quet states form an orthonormal and complete ba- 
sis. 



APPENDIX D: INVERSE OF 

The inverse of (G™)„i„(cj) is calculated for any and 
A{t) as follows. First, using expression ([TT]) we have 

G^" ^ = Afc • Qf,^{uj) ■ since A^ is unitary as proved 
in Appendix O Among the terms in the diagonal matrix 
Q^^{uj), [uj + fi - (efc)o + i'n]Srnn commutcs with A^, giv- 
ing a trivial result. The only nontrivial part, nQSmm is 
evaluated as 



Y^iAkUm {Ai)in=f Y^S^n.-ny)-U(.-y) exp(^-^ f^dz [.fc_,^ 



dx 

-IT J -TT 

dx r 



(z/n) - (efc)o] 



dy in[dj[x - y)]e*(™-»«) exp ( - - y dz [efc_eA(./o) - 



f dx f 

/ dy d{x - y)[mn - ek_eA{x/n) + {'^k)a]e' 

J -TT J -TT 

X exp ^-^ J dz [efe_eA(z/a) - (efe)o]^ = [mfl + {ek)o]Smn - (e/e)™-™- 



(Dl) 



Between the second and the third lines we have integrated 
by parts. The contribution coming from the boundary is 
negligible due to the presence of the delta function. Thus 
the simple expression [Eq. results. 

We note that relation (jPip for the Floquet indices 
TO = n = reduces to the geometrical phase, 7^ = 
I-i'^/2 '"fc(0*^t'"fc(^)- From Eq. (jDl[) . we can show that 
it is exactly absent in the single-band noninteracting sys- 
tem, 7fc = -2TTY,eJ{^k)oi{^l)m = 0. 



APPENDIX E: JOINT DENSITY OF STATES IN 
ARBITRARY FINITE DIMENSIONS 



An analytic expression for JDOS defined by Eq. ([221) 
can be obtained in arbitrary dimensions. We first sub- 
stitute the delta functions in Eq. ((22|) with the integrals 
over auxiliary variables s and s: 



r°° r°° dn 

^ \ ^ ^2ity^^(s cos ki-\-s sink i) 
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FIG. 12: The JDOS p{() for d = 3 and d = 00 in units of f . 



We then replace e and e with C and 9 in accordance with 
e = ^cos9 and e = Csin0, and change the integral vari- 
ables as s = ^ sin and s = ^ cos (f>. After the Integra- 
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tions, we obtain 



P(C) 



27r 



(El) 



This is the general expression for the JDOS in d di- 
mensions. Note that the JDOS is independent of 9 in 
any dimension. The infinite-dimensional case is read- 
ily reproduced since Jo{z) = 1 — (z/2)^ -I- O(z^) and 
t = t*/2Vd, the factor [Jo{2t^)]'^ converges to e'^^/^)" 
in the hmit d — > oo. With a formula for the integral of 
the Besscl function, we reproduce the known JDOS in 
infinite dimensions)^ 



1 



(E2) 



In the case of finite dimensions, we can systematically 
deduce the JDOS from Eq. (|Eip with an appropriate in- 
tegral formula related to the Bessel function. In the fol- 
lowing we list the derived expression of the JDOS for 
d= 1,2,3: 



Pd=2{C) 



27r(2t) 



SiC-2t), 
1 



< C < 4t 
4t < C 



(E3) 



(E4) 



P<i=3(C) 



1 



K 4, 



(C + 2t)3(6t-C) 



27r3/(2iFC 



K 



1 / (C + 2t)3(6t-C) 
4V (2i)3C 







< C < 2t 



2t<C<6t 



6t<C 



(E5) 



where K{k) is the elliptic integral of the first kind. In Fig. [T2]we plot the JDOS for d = 3 and d = oo. We can observe 
that pd=3 diverges at C = 2t, which originates from the van Hove singularity of the simple cubic lattice. 
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